Robust charge-density-wave correlations in the electron-doped single-band Hubbard model

There is growing evidence that the hole-doped single-band Hubbard and t − J models do not have a superconducting ground state reflective of the high-temperature cuprate superconductors but instead have striped spin- and charge-ordered ground states. Nevertheless, it is proposed that these models may still provide an effective low-energy model for electron-doped materials. Here we study the finite temperature spin and charge correlations in the electron-doped Hubbard model using quantum Monte Carlo dynamical cluster approximation calculations and contrast their behavior with those found on the hole-doped side of the phase diagram. We find evidence for a charge modulation with both checkerboard and unidirectional components decoupled from any spin-density modulations. These correlations are inconsistent with a weak-coupling description based on Fermi surface nesting, and their doping dependence agrees qualitatively with resonant inelastic x-ray scattering measurements. Our results provide evidence that the single-band Hubbard model describes the electron-doped cuprates.

Given their ubiquity, these CDW correlations must be accounted for by any proposed effective model for the cuprates. Evidence for charge modulations, both in the form of unidirectional stripe correlations or short-range CDW correlations, has now been found in a variety of finite-temperature quantum Monte-Carlo (QMC) simulations of the h-doped Hubbard model 13,[34][35][36][37][38] . These simulations are generally restricted to high temperatures by the Fermion sign problem 10,34,35 (except for very recent constrained path QMC calculations 38 ) and focus on the h-doped model. However, the observed cuprate CDWs exhibit a significant electron-hole asymmetry. On the h-doped side, they can intertwine with spin-density modulations to form stripes while they coexist with uniform antiferromagnetic (AFM) correlations on the e-doped side 25,28,30 . These differences have raised questions on whether the e-and h-doped CDWs share a common origin 29,30 . Addressing these questions requires detailed calculations for the edoped Hubbard model, which analyze the spin and charge correlations. To date, such calculations have not been performed.
Here, we provide insights into this question by studying the spin and charge correlations of the e-doped two-dimensional Hubbard model and contrasting them with the h-doped case, using the dynamical cluster approximation (DCA) 39 and a nonperturbative QMC cluster solver 40 (see Methods). Working on large (16 × 4) rectangular clusters that are wide enough to support large-period stripe correlations if they are present 34 , we vary the electron density 〈n〉 across both sides of the phase diagram to contrast the correlations in each case. Our calculations uncover robust two-component CDW correlations on the e-doped side, which consists of superimposed checkerboard (0.5, 0.5) (in the unit of 2π/a implied for all vectors in k space) and unidirectional Q CDW = ( ± δ c , 0) components. These CDW correlations appear to be decoupled from any stripe-like modulations of the spins and instead coexist with short-range antiferromagnetic correlations. This behavior is in direct contrast to the h-doped case, where we find evidence for fluctuating stripe correlations in both the charge and spin degrees of freedom 34 . Our results agree with experimental observations on the e-doped cuprates, including the observed doping dependence of Q CDW . This supports the notion that the single-band Hubbard model captures the e-doped side of the high-T c phase diagram. Figure 1 compares the static (ω = 0) charge χ c (Q, 0) and spin χ s (Q, 0) susceptibilities for the h-(〈n〉 = 0.8) and e-doped (〈n〉 = 1.2) Hubbard model with U=t = 6, t 0 =t = À 0:2, and varying temperature (see Supplementary Note 1 for error information of the lowest temperature results). In the h-doped case (Fig. 1a, c, e), unidirectional charge and spin stripes form as the temperature is lowered, consistent with prior finite-temperature studies 34,35,38 . These correlations manifest as incommensurate peaks in the static susceptibility at wave vectors Q c = ( ± δ c , 0) and Q s = (0.5 ± δ s , 0.5) for the charge and spin channels, respectively. For the spin channel in Fig. 1e, the dashed line shows a fit of the lowest temperature data using two Lorentzian functions. Figure 1c plots the charge correlations for the h-doped case along the (Q x , 0.5) direction, where we observe a weak double-peak structure emerging at the lowest accessible temperature. This modulation is weaker than the Q c structure in Fig. 1a, such that the charge correlations are predominantly stripe-like.

Results
We observe qualitatively different correlations in the e-doped case shown in Fig. 1b, d, f. At high temperature (β ≤ 6/t), χ c (Q, 0) has a single broad peak centered at q = (0, 0), which can again be decomposed into two incommensurate Lorentzian peaks centered at ± δ c , indicative of a unidirectional charge stripe. As the temperature is lowered, these peaks sharpen and become discernible without fits while the q-independent background remains constant. The charge correlations also have a relatively temperature-independent (0.5, 0.5) component ( Fig. 1e) of similar strength as the stripe-like charge correlations. In contrast to the h-doped case, we find no indication of a spin-stripe at this doping; the spin susceptibility has a single peak centered at (0.5, 0.5) for all accessible temperatures.
Comparing panels a and b, we see that the charge stripe correlations in the h-and e-doped systems develop differently as the temperature decreases. In the h-doped case (Fig. 1a), the incommensurate peaks grow while δ c shifts to smaller values as the temperature decreases 34 . In the e-doped case, the incommensurate peaks grow (Fig. 1b) while the value of χ c (Q, 0) near zone center is suppressed, resulting in well-defined peaks centered at ≈(±0.25, 0). In addition, the (0.5, 0.5) component is significantly stronger in the e-doped case (Fig. 1d). However, since the height of the incommensurate peaks in Fig. 1b does not appear to level off, the stripe correlations could dominate at lower temperatures.
The corresponding correlation functions in real-space, obtained at the lowest accessible temperatures (βt = 6 and 10 for the h-and edoped cases, respectively) are plotted in  Figure 2c, d show the staggered spin correlation function for the h-and e-doped cases, respectively. Here, the blue region in the middle represents one AFM domain, while the red region on both sides indicates neighboring AFM domains with the opposite phase. (Note that the staggered spin correlation function that is plotted contains an additional negative sign on the B-sublattice as explained in the Methods section.) Consistent with Fig. 1 and as observed before 34 , the hdoped case has a clear stripe pattern whose period is roughly twice that of the charge modulations. In contrast, the e-doped model is dominated by short-range AFM correlations with only a single phase inversion appearing at longer distances. To summarize, the h-doped system has intertwined spin and charge stripe correlations, while the edoped system manifests CDW correlation with stripe-like Q = (δ c , 0) and checkerboard (0.5, 0.5) components and nearly uniform AFM spin correlations. Figure 3 examines how the stripe component of the charge correlations develops for the e-doped case with density and t 0 for a fixed inverse temperature β = 8/t. Figure 3a shows χ c (Q, 0) along the (Q x , 0) direction for various densities, while holding t 0 =t = À 0:2 fixed. At 〈n〉 = 1.125, the charge susceptibility has a single peak centered at Q = 0. With further electron doping, the peak splits into two well-defined peaks that become discernible without fitting for 〈n〉 ≥ 1.2. At the same time, the uniform background grows with doping due to the increased metallicity in the system 35 . Figure 3b presents the same quantity for various values of t 0 , while fixing 〈n〉 = 1.2. At t 0 =t = À 0:1, the charge susceptibility again has a single broad peak at Q = 0. As |t 0 | increases, the middle peak is suppressed, leading to the appearance of a pair of incommensurate peaks. At the same time, the uniform background remains almost unchanged.
To assess these trends quantitatively, we fit the curves in Fig. 3a with three Lorentzian functions centered at Q x = 0 and ± Q c plus a constant background to extract the wave vector δ c = Q c . As shown in Fig. 3c, δ c increases approximately linearly as a function of the doping ρ = 〈n〉 − 1. This trend agrees with experimental observations for the edoped cuprates 28 ; however, the observed Q cdw (ρ) curve is shifted to lower doping levels relative to our results. We also extracted δ c as a function of t 0 , as shown in Fig. 3d, where we find that δ c linearly shifts to larger values with |t 0 |. Extrapolating the experimental data in Fig. 3c to ρ = 0.2 yields Q cdw ≈ 0.32 r.l.u., which corresponds to |t 0 =t| ≈ 0:4 in our model. This value is much larger than the typical values used to model the e-doped cuprates using the single-band model.
RIXS experiments have found that the CDW in the e-doped cuprates has a significant dynamical component 25,28,30 . To compare with these measurements, we show in Fig. 4 the dynamical spin and charge structure functions S(Q, ω) and N(Q, ω) obtained using a parameter-free differential evolution analytic continuation (DEAC) algorithm 41 . (A comparison of the DEAC results to those obtained with more conventional techniques is provided in Supplementary Note 3.) The dynamic spin structure factor S(Q, ω) displays the typical persistant antiferromagnetic paramagnon spectrum obtained with other QMC methods 42 , with large spectral weight at Q = (0.5, 0) at an energy near ω = t and a downward dispersion towards Q = 0. The dynamic charge structure factor exhibits similar behavior with a large spectral weight near the zone edge and a downward dispersion towards the zone center. Still, the spectral weight is concentrated at higher energies.
As an approximation of the predicted RIXS intensity, Fig. 4c plots a superposition of S(Q, ω) and N(Q, ω). Due to the vastly different energy scales in the spin and charge correlations, one sees two distinct upward dispersing branches, i.e., a low-energy branch due to the spin correlations and a high-energy branch due to the charge correlations. Since the spin correlations have a much larger amplitude, the lowenergy branch has a stronger overall intensity, while the higher energy (charge) branch is muted. The behavior for the dynamical correlations near zone center agrees well with the data reported in ref. 30. Our results for the higher energy charge excitations call for future RIXS measurements in this regime; however, the high-energy portion of the charge excitations may overlap with the intra-atomic orbital excitations on the Cu atom (the so-called dd-excitations), which are commonly observed at the Cu L-edge 43 .

Discussion
The correlations in the e-doped case cannot be attributed solely to the weak-coupling Lindhard physics 44 . For example, as discussed in Supplementary Note 4, if we adjust the effective U in a random-phase approximation (RPA) to match the charge correlations χ c (Q c , 0) then the predicted correlations at (0.5, 0.5) are much stronger than those reported here. Similarly, we do not resolve any peak structure in χ s (Q, 0) along the (Q x , 0) direction as one would expect in such a weakcoupling framework. While the peak positions predicted by RPA are close to the values reported here for 〈n〉 = 1.2, the resulting temperature, doping, and t 0 evolution is inconsistent with the observed behavior (see Supplementary Note 4).
We have demonstrated that the CDW correlations observed in DCA simulations of the single-band Hubbard model have a pronounced particle-hole asymmetry. Our results for the e-doped system are also in qualitative agreement with RIXS experiments on Nd 2−x Ce x CuO 4 25,30 ; however, there are some notable quantitative differences. For example, our predicted δ c (ρ) dependence is shifted to higher electron doping in comparison to experiment. This   χ c (r, 0) for the e-doped system (t 0 = À 0:2t,hni = 1:2) at the lowest accessible inverse temperature βt = 10. b χ c (r, 0) for the h-doped system (t 0 = À 0:2t,hni = 0:8) at the lowest accessible inverse temperature βt = 6. c, d show the staggered spin-spin correlations χ s,stag (r, 0) for the h-and e-doped cases, respectively. + and − signs indicate the sign of the correlations whose absolute mean is larger than two standard errors. The dashed lines indicate the approximate nodes in the spin and charge stripe modulations. discrepancy may be related to challenges in determining the carrier concentration in the CuO 2 planes. ARPES measurements of Pr 1.3−x La 0.7 Ce x CuO 4−δ (PLCCO), for example, have suggested that the doped electron concentration of CuO 2 plane can be larger than the Ce concentration x by up to 0.08 e/Cu, depending on the annealing method 45 . This discrepancy is comparable to the shift observed in Fig. 3c. Adjusting the value of t 0 could also partially account for this difference; single-band fits to the measured Fermi surface of PLCCO estimate |t 0 =t|≈0:34 À 0:4. Finally, the periodicity of the charge modulations may be affected by the DCA mean-field 34 . In the h-doped case, DCA and DQMC predict different stripe periods for the same model parameters. Nevertheless, our calculations reproduce the qualitative doping dependence observed in the real materials and support the notion that the single-band Hubbard model captures the physics of the e-doped cuprates.

The model
We consider the single-band Hubbard model on a two-dimensional square lattice. The Hamiltonian is Here, c y i,σ (c i,σ ) creates (annihilates) a spin-σ ( = ↑, ↓) electron at site i = a(i x , i y ), where a = 1 is the lattice constant, t i,j is the hopping integral between sites i and j, μ is the chemical potential, and U is the Hubbard repulsion. To model the cuprates, we set t i,j = t and t 0 for nearest and next-nearest-neighbor hopping, respectively, and zero otherwise, and take U/t = 6 unless stated otherwise.

The dynamical cluster approximation
We study the model in Eq. (1) using DCA 8,39,46 as implemented in the DCA++ code 47 . The DCA represents the infinite bulk lattice in the thermodynamic limit by a finite-size cluster embedded in a selfconsistent dynamical mean-field. The intra-cluster correlations are treated exactly, while the mean-field approximates the inter-cluster degrees of freedom. We use rectangular N c = 16 × 4 clusters that are large enough to support spin and charge stripe correlations if they are present in the model 34 .
Assuming short-ranged correlations, the self-energy Σ(k, iω n ) can be approximated by the cluster self-energy Σ(K, iω n ), where K are the cluster momenta. We obtain the coarse-grained single-particle Green by averaging the lattice Green function G(k, iω n ) over the N/N c momenta k 0 in a patch about the cluster momentum K. (N and N c are the number of site in the lattice and cluster, respectively.) In this way, the bulk problem is reduced to a finite-size-cluster problem. We solve the cluster problem with the continuous-time, auxiliaryfield quantum Monte-Carlo algorithm (CT-AUX) 40,48 . The expansion order of the CT-AUX QMC algorithm is typically 500-1500, depending on the temperature and value of t 0 . Depending on the value of the average fermion sign for a given parameter set, we measure 1−5 × 10 8 samples for the correlation functions. Six to eight iterations of the DCA loop are typically needed to obtain good convergence for the DCA cluster self-energy.
To study the spin and charge correlations, we measure the static (ω = 0) real-space spin-spin correlation function and density-density correlation function Here, r = (r x , r y )a and i = (i x , i y )a are lattice vectors andŜ z i = 1 2 ðn i," À n i,# Þ and n i = P σ c y i,σ c i,σ are the local spin-z and density operators of the effective cluster problem. The correlation functions χ c,s (r, 0) are measured directly by the cluster solver. The static momentum-space susceptibilities χ c (Q, 0) and χ s (Q, 0) are obtained by a Fourier transform. We also plot the staggered spin correlation function χ s,stag ðrÞ = ðÀ1Þ ðr x + r y Þ Sðr, 0Þ to highlight the relative phases of the antiferromagnetic domains.

Analytic continuation
The dynamical charge N(Q, ω) and spin S(Q, ω) structure factors are obtained from the imaginary part of the charge and spin β/t = 10. c shows the sum of the two as a crude approximation for Cu L-edge RIXS spectra. The black circles and diamonds indicate the locations of the maxima in S(Q, ω) and N(Q, ω), respectively. All panels are plotted on the same color scale, as indicated on the right. susceptibilities using NðQ, ωÞ = Imχ c ðQ, ωÞ 1 À e Àβω , ð5Þ and SðQ, ωÞ = Imχ s ðQ, ωÞ 1 À e Àβω : ð6Þ The real frequency susceptibility is related to the imaginary time susceptibility by the integral equation χ s,c ðQ, τÞ = Z 1 0 dω π e Àτω + e ÀðβÀτÞω 1 À e Àβω Im χ s,c ðQ, ωÞ: Inverting this relationship to obtain ImχðQ,ωÞ is an ill-posed problem. We used three independent methods to perform the analytic continuation to gauge our relative confidence in the various features. These include the method of Maximum Entropy 49 , a parameter-free differential evolution algorithm 41 , and stochastic optimization 50 . The results obtained using the differential evolution algorithm are shown in the main text, while the remaining results are provided in Supplementary Note 3.

Data availability
The data in this study are available at https://doi.org/10.5281/zenodo. 7829220.